library_load <- suppressMessages(
list(
# Seurat
library(Seurat),
# Scran doubletFinder implementation
library(scDblFinder),
# FastMNN implementation
library(SeuratWrappers),
# Data
library(tidyverse),
# Reticulate
library(reticulate),
# Biomart
library(biomaRt),
# Plotting
library(ggplot2),
library(knitr),
library(kableExtra),
library(IRdisplay)
)
)
random_seed <- 42
set.seed(random_seed)
options(warn=-1)
# Set working directory to project root
setwd("/research/peer/fdeckert/FD20200109SPLENO")
# Source files
source("plotting_global.R")
source("bin/SeuratQC.R")
# SingleRPlot
source("bin/SingleRQC.R")
# Files
so_raw_file <- "data/object/raw.rds"
so_qc_rds <- "data/object/qc.rds"
so_qc_h5ad <- "data/object/qc.h5ad"
# Plotting Theme
ggplot2::theme_set(theme_global_set()) # From project global source()
so_raw <- readRDS(so_raw_file)
so_raw$tissue <- factor(so_raw$tissue, levels=c("Myeloid", "Progenitor"))
so_raw$treatment <- factor(so_raw$treatment, levels=c("NaCl", "CpG"))
Empty droplets were determined with CellRanger V3.0.2 Lun et al., 2019 EmptyDrop heuristic. RNAse activity of granulocytes might be wrongly identified as empty cells by CellRanger.
Typical Sample A steep drop-off is indicative of good separation between the cell-associated barcodes and the barcodes associated with empty GEMs. A ideal barcode rank plot has a distincitve shape, which is referred to as a "cliff and knee".
Heterogeneous Sample Heterogeneous populations of cells in a sample result in two "cliff and knee" distributions. However, there should still be clear separation between the bacodes.
Compromised Sample Round curve and lack of steep cliff may indicate low sample quality or loss of single-cell behavior. This can be due to a wetting failure, premature cell lysis, or low cell viability.
Compromised Sample Defined cliff and knee, but the total number of barcodes detected may be lower than expected. This can be caused by a sample clog or inaccurate cell count.
options(repr.plot.width=10, repr.plot.height=5)
rank_plot_qc_1 <- rank_plot_qc(so_raw, color_color=color$cellranger_class, formular=tissue~sample_rep+treatment)
rank_plot_qc_1
ggsave(rank_plot_qc_1, filename="result/qc/plot/rank_plot_qc.png", width=10, height=5)
so_qc <- subset(so_raw, subset=cellranger_class == "Cell")
# Get Seurat cell cycle genes
cc_genes_seurat_s <- cc.genes.updated.2019$s.genes
cc_genes_seurat_g2m <- cc.genes.updated.2019$g2m.genes
# Get mouse orthologs from human gene simbols
httr::set_config(httr::config(ssl_verifypeer=FALSE))
hgnc_mart <- useMart("ensembl", dataset="hsapiens_gene_ensembl", host="https://dec2021.archive.ensembl.org/")
mm_mart <- useMart("ensembl", dataset="mmusculus_gene_ensembl", host="https://dec2021.archive.ensembl.org/")
cc_genes_seurat_s <- getLDS(attributes=c("hgnc_symbol"), filters="hgnc_symbol", values=cc_genes_seurat_s, mart=hgnc_mart, attributesL=c("mgi_symbol"), martL=mm_mart, uniqueRows=TRUE)
cc_genes_seurat_s <- cc_genes_seurat_s[, 2]
cc_genes_seurat_g2m <- getLDS(attributes=c("hgnc_symbol"), filters="hgnc_symbol", values=cc_genes_seurat_g2m, mart=hgnc_mart, attributesL=c("mgi_symbol"), martL=mm_mart, uniqueRows=TRUE)
cc_genes_seurat_g2m <- cc_genes_seurat_g2m[, 2]
# Compute cell cycle score
so_qc <- CellCycleScoring(so_qc, s.features=cc_genes_seurat_s, g2m.features=cc_genes_seurat_g2m, set.ident=FALSE)
colnames(so_qc@meta.data) <- gsub("Phase", "cc_phase_class", colnames(so_qc@meta.data))
colnames(so_qc@meta.data) <- gsub("S.Score", "msS_RNA", colnames(so_qc@meta.data))
colnames(so_qc@meta.data) <- gsub("G2M.Score", "msG2M_RNA", colnames(so_qc@meta.data))
so_qc$msCC_diff_RNA <- so_qc$msS_RNA - so_qc$msG2M_RNA
so_qc$cc_phase_class <- factor(so_qc$cc_phase_class, levels=names(color$cc_phase_class))
SingleR identifies marker genes from the reference and uses them to compute assignment score (based on the Spearman correlation across markers) for each cell in the test dataset against each label in the reference. The label with the highest doublet_stat is the assigned to the test cell, possibly with further fine-tuning to resolve closely related labels.
first.labels: Labels before fine-tuning
labels: Labels after fine-tuning
pruning: labels after pruning
if(FALSE) {
# Load reference data
ref <- ImmGenData()
# Seurat object to SingleCellExperiment
sce <- SingleCellExperiment(list(counts=so_qc@assays$RNA@counts))
# # Predict labels
label_main <- SingleR::SingleR(test=sce, ref=ref, labels=ref$label.main, assay.type.test="counts", de.method="classic")
label_fine <- SingleR::SingleR(test=sce, ref=ref, labels=ref$label.fine, assay.type.test="counts", de.method="classic")
saveRDS(label_main, "result/SingleR/label_main_immgen.rds")
saveRDS(label_fine, "result/SingleR/label_fine_immgen.rds")
} else {
label_main <- readRDS("result/SingleR/label_main_immgen.rds")
label_fine <- readRDS("result/SingleR/label_fine_immgen.rds")
}
# Add labels to Seurat object
label_main_meta <- as.data.frame(label_main) %>% dplyr::select(pruned.labels, tuning.scores.first) %>% dplyr::rename(label_main_immgen=pruned.labels, delta_score_main_immgen=tuning.scores.first)
label_fine_meta <- as.data.frame(label_fine) %>% dplyr::select(pruned.labels, tuning.scores.first) %>% dplyr::rename(label_fine_immgen=pruned.labels, delta_score_fine_immgen=tuning.scores.first)
so_qc <- AddMetaData(so_qc, cbind(label_main_meta, label_fine_meta))
if(FALSE) {
# Load reference data
ref <- readRDS("data/haemosphere/se_haemosphere.rds")
# Seurat object to SingleCellExperiment
sce <- SingleCellExperiment(list(counts=so_qc@assays$RNA@counts))
# # Predict labels
label_main <- SingleR::SingleR(test=sce, ref=ref, labels=ref$label.main, assay.type.test="counts", de.method="classic")
label_fine <- SingleR::SingleR(test=sce, ref=ref, labels=ref$label.fine, assay.type.test="counts", de.method="classic")
saveRDS(label_main, "result/SingleR/label_main_haemosphere.rds")
saveRDS(label_fine, "result/SingleR/label_fine_haemosphere.rds")
} else {
label_main <- readRDS("result/SingleR/label_main_haemosphere.rds")
label_fine <- readRDS("result/SingleR/label_fine_haemosphere.rds")
}
# Add labels to Seurat object
label_main_meta <- as.data.frame(label_main) %>% dplyr::select(pruned.labels, tuning.scores.first) %>% dplyr::rename(label_main_haemosphere=pruned.labels, delta_score_main_haemosphere=tuning.scores.first)
label_fine_meta <- as.data.frame(label_fine) %>% dplyr::select(pruned.labels, tuning.scores.first) %>% dplyr::rename(label_fine_haemosphere=pruned.labels, delta_score_fine_haemosphere=tuning.scores.first)
so_qc <- AddMetaData(so_qc, cbind(label_main_meta, label_fine_meta))
options(repr.plot.width=10, repr.plot.height=10)
ggplot(so_qc@meta.data, aes(y=label_fine_haemosphere, fill=label_fine_haemosphere)) +
geom_bar() +
xlab("log10(cell count)") + ylab("") + ggtitle("Fine Haemosphere label") +
scale_x_log10(breaks=c(0, 10, 100, 1000)) +
facet_grid(tissue~treatment+sample_rep) +
scale_fill_manual(values=color$label_fine_haemosphere)
We recognize a substantial number or cells annotated as erythrocytes in the Myeloid sort. Since we do not expect those cells based on the FACS panel they are probably doublet cells or phagocytes which took up a erythrocyte. As shown in the scatter plots below those cells show either a very high or low feature count respectively. While it might be of interest keeping the phagocytes the low number of features detected – probably to the exaggerated Hb mRNA count – makes downstream analysis difficult. We therefore exclude cells from the Myeloid sort which were annotated as Ery.
qc_class_set <- function(so) {
# QC filter for all cells
so$pMt_RNA_min <- -Inf
so$pMt_RNA_max <- 5
so$nCount_RNA_min <- 1500
so$nCount_RNA_max <- max(so$nCount_RNA)
so$qc_class <- ifelse(
so$cellranger_class == "Cell" &
so$nCount_RNA <= so$nCount_RNA_max &
so$nCount_RNA > so$nCount_RNA_min &
so$pMt_RNA <= so$pMt_RNA_max &
so$pMt_RNA > so$pMt_RNA_min,
"pass", "fail"
)
# Remove Ery from Myeloid sort before computing nFeature_RNA
if(so$tissue[1]=="Myeloid") {
so$qc_class <- ifelse(
so$label_main_haemosphere=="Ery" | so$label_fine_haemosphere %in% c("preCFUE", "CFUE", "pbEry", "poEry", "Retic") | so$pHb_RNA >= 5,
"fail", so$qc_class
)
}
so_tmp <- subset(so, subset=qc_class=="pass")
so$nFeature_RNA_max <- quantile(so_tmp$nFeature_RNA, 0.99)
so$nFeature_RNA_min <- quantile(so_tmp$nFeature_RNA, 0.01)
so$qc_class <- ifelse(
so$qc_class == "pass" &
so$nFeature_RNA <= so$nFeature_RNA_max &
so$nFeature_RNA > so$nFeature_RNA_min,
"pass", "fail"
)
return(so)
}
so_qc <- Seurat::SplitObject(so_qc, split.by="sample_name")
so_qc <- lapply(so_qc, qc_class_set)
so_qc <- merge(so_qc[[1]], so_qc[2:length(so_qc)])
so_qc$tissue <- factor(so_qc$tissue, levels=c("Myeloid", "Progenitor"))
so_qc$treatment <- factor(so_qc$treatment, levels=c("NaCl", "CpG"))
density_plot_qc_1 <- density_plot_qc(so=so_qc, title="Density plot UMI count", x=nCount_RNA, xlab="log10(UMI count)", min=nCount_RNA_min, max=nCount_RNA_max, fill_color=color$tissue)
density_plot_qc_2 <- density_plot_qc(so=so_qc, title="Density plot Feature count", x=nFeature_RNA, xlab="log10(Feature count)", min=nFeature_RNA_min, max=nFeature_RNA_max, fill_color=color$tissue)
density_plot_qc_3 <- density_plot_qc(so=so_qc, title="Density plot Mt %", x=pMt_RNA, xlab="Mt [%]", min=0, max=pMt_RNA_max, log10=FALSE, fill_color=color$tissue)
ggsave(density_plot_qc_1, filename="result/qc/plot/density_plot_qc_umi.png", width=10, height=5)
ggsave(density_plot_qc_2, filename="result/qc/plot/density_plot_qc_feature.png", width=10, height=5)
ggsave(density_plot_qc_3, filename="result/qc/plot/density_plot_qc_mt.png", width=10, height=5)
options(repr.plot.width=30, repr.plot.height=5)
density_plot_qc_1 + density_plot_qc_2 + density_plot_qc_3 + plot_layout(nrow=1) & theme(legend.position="bottom")
scattern_plot_qc_1 <- scattern_plot_qc(so=so_qc, title="Mitochondrial gene percentage", color=pMt_RNA)
scattern_plot_qc_2 <- scattern_plot_qc(so=so_qc, title="Hemoglobin gene percentage", color=pHb_RNA)
scattern_plot_qc_3 <- scattern_plot_qc(so=so_qc, title="Ribsonmal gene percentage", color=pRb_RNA)
scattern_plot_qc_4 <- scattern_plot_qc(so=so_qc, title="Fine labels", color=label_fine_haemosphere) + scale_color_manual(values=color$label_fine_haemosphere)
scattern_plot_qc_5 <- scattern_plot_qc(so=so_qc, title="QC class", color=qc_class) + scale_color_manual(values=c("fail"="red", "pass"="black"))
ggsave(scattern_plot_qc_1, filename="result/qc/plot/sc_mt.png", width=10, height=5)
ggsave(scattern_plot_qc_2, filename="result/qc/plot/sc_hg.png", width=10, height=5)
ggsave(scattern_plot_qc_3, filename="result/qc/plot/sc_rb.png", width=10, height=5)
ggsave(scattern_plot_qc_4, filename="result/qc/plot/sc_fl.png", width=10, height=5)
ggsave(scattern_plot_qc_5, filename="result/qc/plot/sc_qc.png", width=10, height=5)
options(repr.plot.width=15, repr.plot.height=15)
scattern_plot_qc_1 + scattern_plot_qc_2 + scattern_plot_qc_3 + scattern_plot_qc_4 + scattern_plot_qc_5 + plot_layout(nrow=3) & theme(legend.position="bottom")
box_plot_qc_1 <- box_plot_qc(so=so_qc, y=nCount_RNA, fill=tissue, ylab="UMI [count]", ymin=0)
box_plot_qc_2 <- box_plot_qc(so=so_qc, y=nFeature_RNA, fill=tissue, ylab="Feature [count]", ymin=0)
box_plot_qc_3 <- box_plot_qc(so=so_qc, y=pMt_RNA, fill=tissue, ylab="Mt [%]", ymin=0, ymax=100)
box_plot_qc_4 <- box_plot_qc(so=so_qc, y=pHb_RNA, fill=tissue, ylab="Hb [%]", ymin=0, ymax=100)
box_plot_qc_5 <- box_plot_qc(so=so_qc, y=pRb_RNA, fill=tissue, ylab="Rbl [%]", ymin=0, ymax=100)
ggsave(box_plot_qc_1[[1]] + box_plot_qc_1[[2]] + plot_layout(guides="collect", ncol=1) & theme(legend.position="bottom"), filename="result/qc/plot/box_plot_qc_1.png", width=7.5, height=5)
ggsave(box_plot_qc_2[[1]] + box_plot_qc_2[[2]] + plot_layout(guides="collect", ncol=1) & theme(legend.position="bottom"), filename="result/qc/plot/box_plot_qc_2.png", width=7.5, height=5)
ggsave(box_plot_qc_3[[1]] + box_plot_qc_3[[2]] + plot_layout(guides="collect", ncol=1) & theme(legend.position="bottom"), filename="result/qc/plot/box_plot_qc_3.png", width=7.5, height=5)
ggsave(box_plot_qc_4[[1]] + box_plot_qc_4[[2]] + plot_layout(guides="collect", ncol=1) & theme(legend.position="bottom"), filename="result/qc/plot/box_plot_qc_4.png", width=7.5, height=5)
ggsave(box_plot_qc_5[[1]] + box_plot_qc_5[[2]] + plot_layout(guides="collect", ncol=1) & theme(legend.position="bottom"), filename="result/qc/plot/box_plot_qc_5.png", width=7.5, height=5)
options(repr.plot.width=15, repr.plot.height=20)
box_plot_qc_1[[1]] + box_plot_qc_1[[2]] + box_plot_qc_2[[1]] + box_plot_qc_2[[2]] +
box_plot_qc_3[[1]] + box_plot_qc_3[[2]] + box_plot_qc_4[[1]] + box_plot_qc_4[[2]] +
box_plot_qc_5[[1]] + box_plot_qc_5[[2]] + plot_spacer() + plot_spacer() + plot_layout(guides="collect", ncol=2) & theme(legend.position="none")
so_qc <- subset(so_qc, subset = qc_class=="pass")
so_qc <- NormalizeData(so_qc, assay="RNA")
so_qc <- ScaleData(so_qc, assay="RNA", verbose=FALSE)
so_qc <- FindVariableFeatures(so_qc, nfeatures=3000, assay="RNA", verbose=FALSE)
so_qc <- RunPCA(so_qc, npcs=50, assay="RNA", features=NULL, verbose=FALSE)
so_qc <- FindNeighbors(so_qc, dims=1:10, k.param=20, verbose=FALSE)
so_qc <- FindClusters(so_qc, verbose=FALSE, resolution=1, algorithm=1, group.singletons=TRUE)
so_qc <- RunUMAP(so_qc, dims=1:50, n.neighbors=30, verbose=FALSE)
options(repr.plot.width=22.5, repr.plot.height=15)
dplot_1 <- dplot(so_qc, reduction="umap", group_by="seurat_clusters", label=TRUE) + ggtitle("Louvain cluster")
dplot_2 <- dplot(so_qc, reduction="umap", group_by="tissue") + scale_color_manual(values=color$tissue) + ggtitle("Tissue")
dplot_3 <- dplot(so_qc, reduction="umap", group_by="treatment") + scale_color_manual(values=color$treatment) + ggtitle("Treatment")
dplot_4 <- dplot(so_qc, reduction="umap", group_by="cc_phase_class") + scale_color_manual(values=color$cc_phase_class) + ggtitle("Cell cycle")
dplot_5 <- dplot(so_qc, reduction="umap", group_by="label_main_haemosphere") + scale_color_manual(values=color$label_main_haemosphere) + ggtitle("Haemosphere label (main)")
dplot_6 <- dplot(so_qc, reduction="umap", group_by="label_fine_haemosphere") + scale_color_manual(values=color$label_fine_haemosphere) + ggtitle("Haemosphere label (fine)")
dplot_7 <- dplot(so_qc, reduction="umap", group_by="sample_rep") + ggtitle("Replicate")
fplot_1 <- fplot(so_qc, reduction="umap", features="Trac") + scale_colour_viridis_c(option="inferno") + ggtitle("Trac")
fplot_2 <- fplot(so_qc, reduction="umap", features="Igkc") + scale_colour_viridis_c(option="inferno") + ggtitle("Igkc")
dplot_comb <- dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + dplot_6 + dplot_7 + fplot_1 + fplot_2 + plot_layout(ncol=3)
dplot_comb
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
so_qc@meta.data %>% dplyr::group_by(seurat_clusters, treatment, sample_rep) %>%
dplyr::summarise(n=n()) %>% data.frame() %>%
tidyr::spread(seurat_clusters, n) %>%
kableExtra::kable("html") %>% as.character() %>% IRdisplay::display_html()
`summarise()` has grouped output by 'seurat_clusters', 'treatment'. You can override using the `.groups` argument.
| treatment | sample_rep | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20 | 21 | 22 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| NaCl | Rep1 | 482 | 22 | 16 | 102 | 120 | 288 | 156 | 112 | 503 | 56 | 300 | 113 | 11 | 128 | 160 | 79 | 67 | 92 | 139 | 42 | 112 | 22 | 8 |
| NaCl | Rep2 | 280 | 82 | 119 | 40 | 113 | 554 | 737 | 565 | 43 | 82 | 170 | 262 | 101 | 135 | 470 | 120 | 229 | 65 | 65 | 62 | 27 | 67 | 29 |
| CpG | Rep1 | 531 | 529 | 239 | 794 | 426 | 131 | 79 | 100 | 457 | 358 | 189 | 164 | 154 | 209 | 38 | 195 | 108 | 78 | 129 | 106 | 90 | 13 | 20 |
| CpG | Rep2 | 154 | 744 | 950 | 310 | 549 | 142 | 127 | 256 | 22 | 483 | 317 | 380 | 569 | 316 | 74 | 316 | 271 | 194 | 64 | 80 | 15 | 1 | 42 |
I use the strategy from Sala et al., 2019 which was also used in Barile et al., 2021 with the source documented here. We only do the refined method with combined samples. The functions used in the script are documented here. We do doublet detection on sample groups (Myeloid and Progenitor cells of a replicated) and integrated those groups per batch and then by treatment with FastMNN
# Split data by sample group
so_qc <- SplitObject(so_qc, split.by="sample_group")
# Compute doublet density per sample group
so_qc <- lapply(so_qc, function(so) {
DefaultAssay(so) <- "RNA"
# Remove empty genes from split
cnt <- GetAssayData(so, assay="RNA", slot="counts")
cnt <- cnt[rowSums(cnt) > 0, ]
# Highly variable features
so <- ScaleData(so, features=rownames(cnt), assay="RNA", verbose=FALSE)
so <- FindVariableFeatures(so, assay="RNA", nfeatures=2000)
# Compute doublet density
so$doublet_score <- scDblFinder::computeDoubletDensity(cnt, subset.row=VariableFeatures(so), size.factors.norm=so$nCount_RNA, dims=50)
so$doublet_score_log2 <- log2(so$doublet_score)
return(so)
}
)
# Compute sub-cluster per sample group
so_qc <- lapply(so_qc, function(so) {
# Set default assay
DefaultAssay(so) <- "RNA"
# Highly variable genes across all data sets
cnt <- GetAssayData(so, assay="RNA", slot="data")
cnt <- cnt[rowSums(cnt)>0, ]
# Cluster function
so <- ScaleData(so, features=rownames(cnt), assay="RNA", verbose=FALSE)
so <- RunPCA(so, npcs=50, features=rownames(cnt), assay="RNA", verbose=FALSE)
so <- FindNeighbors(so, dims=1:50, k.param=15, reduction="pca", verbose=FALSE)
so <- FindClusters(so, verbose=FALSE, resolution=1, algorithm=1)
# Compute sub-cluster
cluster_list <- list()
for(cluster_id in unique(so$seurat_clusters)) {
sink("NUL")
so <- FindSubCluster(so, cluster_id, graph.name="RNA_snn", subcluster.name="doublet_cluster", resolution=1, algorithm=1)
sink()
cluster <- so@meta.data[so$seurat_clusters==cluster_id, ]
cluster <- cluster[, c("cell_id", "doublet_cluster")]
cluster_list[[cluster_id]] <- cluster
}
cluster <- do.call(rbind, cluster_list)
rownames(cluster) <- cluster$cell_id
cluster <- cluster[, "doublet_cluster", drop=FALSE]
so <- AddMetaData(so, cluster)
return(so)
}
)
doublet_stat <- lapply(so_qc, function(so) {
meta <- so@meta.data
# Function to computes MAD for values of x greater the median
mad_upper <- function(x){
x <- x-median(x)
x_mad <- mad(x[x>0], center=0)
return(x_mad)
}
# Comput median doublet density per cluster
doublet_stat <- dplyr::group_by(meta, doublet_cluster, sample_group) %>%
dplyr::summarise(doublet_cluster_density=median(doublet_score))
# Compute median-centered, MAD-estimated variance normal distribution of cluster medians.
doublet_stat$p=pnorm(
doublet_stat$doublet_cluster_density,
mean=median(doublet_stat$doublet_cluster_density),
sd=mad_upper(doublet_stat$doublet_cluster_density),
lower.tail=FALSE
)
return(doublet_stat)
}
)
doublet_stat <- do.call(rbind, doublet_stat)
doublet_stat$fdr=p.adjust(doublet_stat$p, method="fdr")
# Call doublet cluster in meta data
so_qc <- lapply(names(so_qc), function(sample_group){
so <- so_qc[[sample_group]]
so$doublet_class <- ifelse(so$doublet_cluster %in% doublet_stat[doublet_stat$fdr < 0.1 & doublet_stat$sample_group==sample_group, ]$doublet_cluster, "Doublet", "Singlet")
return(so)
}
)
names(so_qc) <- lapply(so_qc, function(so) {so$sample_group[1]})
`summarise()` has grouped output by 'doublet_cluster'. You can override using the `.groups` argument. `summarise()` has grouped output by 'doublet_cluster'. You can override using the `.groups` argument. `summarise()` has grouped output by 'doublet_cluster'. You can override using the `.groups` argument. `summarise()` has grouped output by 'doublet_cluster'. You can override using the `.groups` argument.
# FastMNN on batch
so_qc_1 <- RunFastMNN(so_qc[c("NaCl_Rep2", "NaCl_Rep1")], assay="RNA", reconstructed.assay="MNN")
so_qc_2 <- RunFastMNN(so_qc[c("CpG_Rep2", "CpG_Rep1")], assay="RNA", reconstructed.assay="MNN")
# FastMNN on treatment groups
so_qc <- RunFastMNN(list(so_qc_1, so_qc_2), assay="MNN", reconstructed.assay="MNN")
DefaultAssay(so_qc) <- "MNN"
so_qc <- FindNeighbors(so_qc, dims=1:50, k.param=15, reduction="mnn", graph.name=c("MNN_nn", "MNN_snn"), verbose=FALSE)
so_qc <- FindClusters(so_qc, verbose=FALSE, resolution=1, algorithm=1, graph.name="MNN_snn")
so_qc <- RunUMAP(so_qc, dims=1:50, n.neighbors=30, reduction="mnn", verbose=FALSE)
so_qc$doublet_cluster_int <- so_qc$seurat_clusters
DefaultAssay(so_qc) <- "RNA"
Computing 2000 integration features Computing 2000 integration features Computing 2000 integration features
# Function to computes MAD for values of x greater the median
mad_upper <- function(x){
x <- x-median(x)
x_mad <- mad(x[x>0], center=0)
return(x_mad)
}
doublet_stat <- dplyr::group_by(so_qc@meta.data, doublet_cluster_int) %>%
dplyr::mutate(cluster_count=n()) %>% dplyr::ungroup() %>%
dplyr::group_by(doublet_cluster_int, doublet_class, cluster_count) %>%
dplyr::summarise(doublet_ratio=n()) %>% dplyr::ungroup() %>%
dplyr::mutate(doublet_ratio=doublet_ratio/cluster_count) %>%
dplyr::filter(doublet_class=="Doublet")
doublet_stat$p=pnorm(doublet_stat$doublet_ratio, mean=median(doublet_stat$doublet_ratio), sd=mad_upper(doublet_stat$doublet_ratio), lower.tail=FALSE)
doublet_stat$fdr=p.adjust(doublet_stat$p, method="fdr")
so_qc$doublet_class_int <- ifelse(so_qc$doublet_cluster_int %in% doublet_stat[doublet_stat$fdr < 0.1, ]$doublet_cluster_int, "Doublet", "Singlet")
`summarise()` has grouped output by 'doublet_cluster_int', 'doublet_class'. You can override using the `.groups` argument.
options(repr.plot.width=22.5, repr.plot.height=20)
dplot_1 <- dplot(so_qc, reduction="umap", group_by="treatment") + scale_color_manual(values=color$treatment) + ggtitle("Treatment")
dplot_2 <- dplot(so_qc, reduction="umap", group_by="tissue") + scale_color_manual(values=color$tissue) + ggtitle("Tissue")
dplot_3 <- dplot(so_qc, reduction="umap", group_by="cc_phase_class") + scale_color_manual(values=color$cc_phase_class) + ggtitle("Cell cycle")
dplot_4 <- dplot(so_qc, reduction="umap", group_by="label_fine_haemosphere") + scale_color_manual(values=color$label_fine_haemosphere) + ggtitle("Haemosphere label")
fplot_1 <- fplot(so_qc, reduction="umap", features="doublet_score_log2") + scale_colour_viridis_c(option="inferno") + ggtitle("Log2(Doublet score)")
dplot_5 <- dplot(so_qc, reduction="umap", group_by="doublet_class") + scale_color_manual(values=c("Doublet"="black", "Singlet"="light grey")) + ggtitle("Doublet class sample")
dplot_6 <- dplot(so_qc, reduction="umap", group_by="doublet_class_int") + scale_color_manual(values=c("Doublet"="black", "Singlet"="light grey")) + ggtitle("Doublet class integrated")
fplot_2 <- fplot(so_qc, reduction="umap", features="pMt_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Mt")
fplot_3 <- fplot(so_qc, reduction="umap", features="pHb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Hb")
fplot_4 <- fplot(so_qc, reduction="umap", features="pRb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Rb")
fplot_5 <- fplot(so_qc, reduction="umap", features="nCount_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("UMI count")
fplot_6 <- fplot(so_qc, reduction="umap", features="nFeature_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Feature count")
dplot_comb <- dplot_1 + dplot_2 + dplot_3 + dplot_4 + fplot_1 + dplot_5 + dplot_6 + fplot_2 + fplot_3 + fplot_4 + fplot_5 + fplot_6 + plot_layout(ncol=3)
dplot_comb
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
options(repr.plot.width=22.5, repr.plot.height=5)
dplot_1 <- dplot(so_qc, reduction="umap", group_by="seurat_clusters", label=TRUE) + ggtitle("Louvain cluster")
fplot_1 <- fplot(so_qc, reduction="umap", features="Trac") + scale_colour_viridis_c(option="inferno") + ggtitle("Trac")
fplot_2 <- fplot(so_qc, reduction="umap", features="Igkc") + scale_colour_viridis_c(option="inferno") + ggtitle("Igkc")
dplot_comb <- dplot_1 + fplot_1 + fplot_2 + plot_layout(ncol=3)
dplot_comb
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
so_qc <- subset(so_qc, subset=doublet_class=="Singlet")
so_qc <- subset(so_qc, subset=doublet_class_int=="Singlet")
so_qc <- subset(so_qc, subset=label_fine_haemosphere!="CLP") # Filter out T lymphocytes
so_qc <- subset(so_qc, subset=label_fine_haemosphere!="T-cell") # Filter out T lymphocytes
so_qc <- subset(so_qc, subset=label_fine_haemosphere!="NK") # Filter out T lymphocytes
so_qc <- subset(so_qc, subset=label_main_haemosphere!="NK") # Filter out T lymphocytes
so_qc <- subset(so_qc, subset=seurat_clusters!=18) # Filter out T lymphocytes
so_qc <- subset(so_qc, subset=seurat_clusters!=19) # Filter out B lymphocytes
options(repr.plot.width=22.5, repr.plot.height=20)
dplot_1 <- dplot(so_qc, reduction="umap", group_by="treatment") + scale_color_manual(values=color$treatment) + ggtitle("Treatment")
dplot_2 <- dplot(so_qc, reduction="umap", group_by="tissue") + scale_color_manual(values=color$tissue) + ggtitle("Tissue")
dplot_3 <- dplot(so_qc, reduction="umap", group_by="cc_phase_class") + scale_color_manual(values=color$cc_phase_class) + ggtitle("Cell cycle")
dplot_4 <- dplot(so_qc, reduction="umap", group_by="label_fine_haemosphere") + scale_color_manual(values=color$label_fine_haemosphere) + ggtitle("Haemosphere label")
fplot_1 <- fplot(so_qc, reduction="umap", features="doublet_score_log2") + scale_colour_viridis_c(option="inferno") + ggtitle("Log2(Doublet score)")
dplot_5 <- dplot(so_qc, reduction="umap", group_by="doublet_class") + scale_color_manual(values=c("Doublet"="black", "Singlet"="light grey")) + ggtitle("Doublet class sample")
dplot_6 <- dplot(so_qc, reduction="umap", group_by="doublet_class_int") + scale_color_manual(values=c("Doublet"="black", "Singlet"="light grey")) + ggtitle("Doublet class integrated")
fplot_2 <- fplot(so_qc, reduction="umap", features="pMt_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Mt")
fplot_3 <- fplot(so_qc, reduction="umap", features="pHb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Hb")
fplot_4 <- fplot(so_qc, reduction="umap", features="pRb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Rb")
fplot_5 <- fplot(so_qc, reduction="umap", features="nCount_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("UMI count")
fplot_6 <- fplot(so_qc, reduction="umap", features="nFeature_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Feature count")
dplot_comb <- dplot_1 + dplot_2 + dplot_3 + dplot_4 + fplot_1 + dplot_5 + dplot_6 + fplot_2 + fplot_3 + fplot_4 + fplot_5 + fplot_6 + plot_layout(ncol=3)
dplot_comb
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
options(repr.plot.width=22.5, repr.plot.height=5)
dplot_1 <- dplot(so_qc, reduction="umap", group_by="seurat_clusters", label=TRUE) + ggtitle("Louvain cluster")
fplot_1 <- fplot(so_qc, reduction="umap", features="Trac") + scale_colour_viridis_c(option="inferno") + ggtitle("Trac")
fplot_2 <- fplot(so_qc, reduction="umap", features="Igkc") + scale_colour_viridis_c(option="inferno") + ggtitle("Igkc")
dplot_comb <- dplot_1 + fplot_1 + fplot_2 + plot_layout(ncol=3)
dplot_comb
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
unique(so_qc$label_fine_haemosphere)
unique(so_qc$label_main_haemosphere)
# Store data as RDS
saveRDS(so_qc, so_qc_rds)
# Store data as h5ad
adata <- import("anndata", as="adata", convert=FALSE)
pd <- import("pandas", as="pd", convert=FALSE)
np <- import("numpy", as="np", convert=FALSE)
# Transform dgCMatrix to
X <- GetAssayData(so_qc, assay="RNA", slot="counts") %>% as.matrix() %>% t()
X <- np$array(X, dtype=np$int32)
adata <- adata$AnnData(X=X, obs=so_qc@meta.data)
adata$var_names <- rownames(GetAssayData(so_qc, assay="RNA", slot="counts"))
adata$raw <- adata
adata$write_h5ad(so_qc_h5ad)
None
sessionInfo()
R version 4.1.0 (2021-05-18) Platform: x86_64-conda-linux-gnu (64-bit) Running under: Red Hat Enterprise Linux 8.2 (Ootpa) Matrix products: default BLAS/LAPACK: /home/fdeckert/bin/miniconda3/envs/r.4.1.0-FD20200109SPLENO/lib/libopenblasp-r0.3.15.so locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=en_US.UTF-8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] stats4 grid stats graphics grDevices utils datasets [8] methods base other attached packages: [1] viridis_0.6.2 viridisLite_0.4.0 [3] SingleR_1.8.1 SingleCellExperiment_1.16.0 [5] SummarizedExperiment_1.24.0 Biobase_2.54.0 [7] GenomicRanges_1.46.1 GenomeInfoDb_1.30.1 [9] IRanges_2.28.0 S4Vectors_0.32.4 [11] BiocGenerics_0.40.0 MatrixGenerics_1.6.0 [13] matrixStats_0.62.0 ComplexHeatmap_2.10.0 [15] patchwork_1.1.1 RColorBrewer_1.1-3 [17] IRdisplay_1.1 kableExtra_1.3.4 [19] knitr_1.39 biomaRt_2.50.3 [21] reticulate_1.24 forcats_0.5.1 [23] stringr_1.4.0 dplyr_1.0.9 [25] purrr_0.3.4 readr_2.1.2 [27] tidyr_1.2.0 tibble_3.1.7 [29] ggplot2_3.3.6 tidyverse_1.3.1 [31] SeuratWrappers_0.3.0 scDblFinder_1.8.0 [33] sp_1.4-7 SeuratObject_4.1.0 [35] Seurat_4.1.1 loaded via a namespace (and not attached): [1] rappdirs_0.3.3 pbdZMQ_0.3-7 [3] scattermore_0.8 R.methodsS3_1.8.1 [5] bit64_4.0.5 irlba_2.3.5 [7] DelayedArray_0.20.0 R.utils_2.11.0 [9] data.table_1.14.2 rpart_4.1.16 [11] doParallel_1.0.17 KEGGREST_1.34.0 [13] RCurl_1.98-1.6 generics_0.1.2 [15] ScaledMatrix_1.2.0 cowplot_1.1.1 [17] RSQLite_2.2.13 RANN_2.6.1 [19] future_1.25.0 bit_4.0.4 [21] tzdb_0.3.0 spatstat.data_2.2-0 [23] webshot_0.5.3 xml2_1.3.3 [25] lubridate_1.8.0 httpuv_1.6.5 [27] assertthat_0.2.1 xfun_0.30 [29] hms_1.1.1 evaluate_0.15 [31] promises_1.2.0.1 fansi_1.0.3 [33] progress_1.2.2 dbplyr_2.1.1 [35] readxl_1.4.0 igraph_1.3.1 [37] DBI_1.1.2 htmlwidgets_1.5.4 [39] spatstat.geom_2.4-0 ellipsis_0.3.2 [41] RSpectra_0.16-1 backports_1.4.1 [43] deldir_1.0-6 sparseMatrixStats_1.6.0 [45] vctrs_0.4.1 here_1.0.1 [47] Cairo_1.5-15 remotes_2.4.2 [49] ROCR_1.0-11 abind_1.4-5 [51] batchelor_1.10.0 cachem_1.0.6 [53] withr_2.5.0 progressr_0.10.0 [55] sctransform_0.3.3 prettyunits_1.1.1 [57] scran_1.22.1 goftest_1.2-3 [59] svglite_2.1.0 cluster_2.1.3 [61] lazyeval_0.2.2 crayon_1.5.1 [63] labeling_0.4.2 edgeR_3.36.0 [65] pkgconfig_2.0.3 nlme_3.1-157 [67] vipor_0.4.5 rlang_1.0.2 [69] globals_0.14.0 lifecycle_1.0.1 [71] miniUI_0.1.1.1 filelock_1.0.2 [73] BiocFileCache_2.2.1 modelr_0.1.8 [75] rsvd_1.0.5 rprojroot_2.0.3 [77] cellranger_1.1.0 polyclip_1.10-0 [79] lmtest_0.9-40 Matrix_1.4-1 [81] IRkernel_1.3 zoo_1.8-10 [83] reprex_2.0.1 base64enc_0.1-3 [85] beeswarm_0.4.0 GlobalOptions_0.1.2 [87] ggridges_0.5.3 rjson_0.2.21 [89] png_0.1-7 bitops_1.0-7 [91] R.oo_1.24.0 KernSmooth_2.23-20 [93] Biostrings_2.62.0 blob_1.2.3 [95] DelayedMatrixStats_1.16.0 shape_1.4.6 [97] parallelly_1.31.1 spatstat.random_2.2-0 [99] beachmat_2.10.0 scales_1.2.0 [101] memoise_2.0.1 magrittr_2.0.3 [103] plyr_1.8.7 ica_1.0-2 [105] zlibbioc_1.40.0 compiler_4.1.0 [107] dqrng_0.3.0 clue_0.3-60 [109] fitdistrplus_1.1-8 cli_3.3.0 [111] XVector_0.34.0 listenv_0.8.0 [113] pbapply_1.5-0 MASS_7.3-57 [115] mgcv_1.8-40 tidyselect_1.1.2 [117] stringi_1.7.6 highr_0.9 [119] BiocSingular_1.10.0 locfit_1.5-9.5 [121] ggrepel_0.9.1 tools_4.1.0 [123] future.apply_1.9.0 parallel_4.1.0 [125] circlize_0.4.14 rstudioapi_0.13 [127] uuid_1.1-0 foreach_1.5.2 [129] bluster_1.4.0 metapod_1.2.0 [131] gridExtra_2.3 farver_2.1.0 [133] Rtsne_0.16 digest_0.6.29 [135] BiocManager_1.30.17 rgeos_0.5-9 [137] shiny_1.7.1 Rcpp_1.0.8.3 [139] broom_0.8.0 scuttle_1.4.0 [141] later_1.3.0 RcppAnnoy_0.0.19 [143] httr_1.4.3 AnnotationDbi_1.56.2 [145] colorspace_2.0-3 rvest_1.0.2 [147] XML_3.99-0.9 fs_1.5.2 [149] tensor_1.5 splines_4.1.0 [151] uwot_0.1.11 statmod_1.4.36 [153] spatstat.utils_2.3-0 scater_1.22.0 [155] xgboost_1.6.0.1 plotly_4.10.0 [157] systemfonts_1.0.4 xtable_1.8-4 [159] jsonlite_1.8.0 R6_2.5.1 [161] pillar_1.7.0 htmltools_0.5.2 [163] mime_0.12 glue_1.6.2 [165] fastmap_1.1.0 BiocParallel_1.28.3 [167] BiocNeighbors_1.12.0 codetools_0.2-18 [169] utf8_1.2.2 ResidualMatrix_1.4.0 [171] lattice_0.20-45 spatstat.sparse_2.1-1 [173] curl_4.3.2 ggbeeswarm_0.6.0 [175] leiden_0.3.10 survival_3.3-1 [177] limma_3.50.3 rmarkdown_2.14 [179] repr_1.1.4 munsell_0.5.0 [181] GetoptLong_1.0.5 GenomeInfoDbData_1.2.7 [183] iterators_1.0.14 haven_2.5.0 [185] reshape2_1.4.4 gtable_0.3.0 [187] spatstat.core_2.4-2